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The classical cubic dimer model has a columnar ordering transition that is continuous and described by a crit¬ 
ical Anderson-Higgs theory containing an SU(2)-symmetric complex field minimally coupled to a noncompact 
U(l) gauge theory. Defects in the dimer constraints correspond to monopoles of the gauge theory, with charge 
determined by the deviation from unity of the dimer occupancy. By introducing such defects into Monte Carlo 
simulations of the dimer model at its critical point, we determine the scaling dimensions y 2 = 1.48 ± 0.07 and 
y 2 = 0.20 ± 0.03 for the operators corresponding to defects of charge q = 2 and 3 respectively. These results, 
which constitute the first direct determination of the scaling dimensions, shed light on the deconfined critical 
point of spin-4 quantum antiferromagnets, thought to belong to the same universality class. In particular, the 
positive value of y 3 implies that the transition in the JQ model on the honeycomb lattice is of first order. 


I. INTRODUCTION 

One of the most remarkable consequences of the theory of 
critical phenomena is universality, the occurrence of identi¬ 
cal values for nontrivial measurable quantities at phase tran¬ 
sitions in quite different physical systems (DEI- For Lan¬ 
dau transitions, whose critical properties are described by the 
long-wavelength fluctuations of an order parameter, two tran¬ 
sitions usually belong in the same universality class when 
their space(-time) dimensionalities are the same, and when 
their order parameters have identical symmetry properties at 
the critical point. 

In certain unconventional phase transitions, the order pa¬ 
rameter is not the basic field describing the critical properties, 
but can instead be expressed as a compound object in terms of 
the basic fields 0. Universality can be particularly striking in 
such cases, uniting transitions where the microscopic models 
and phenomenology bear little resemblance. 

The focus of the current work is the noncompact CP ] 
(NCCP 1 ) universality class, proposed as describing a number 
of unconventional phase transitions. These include: the Neel- 
valence-bond solid (VBS) transition 0 in the JQ model 0 
(a quantum antiferromagnet with frustration favoring singlet 
dimers) on certain two-dimensional (2D) lattices; an order¬ 
ing transition in the 3D Heisenberg model with suppression 
of “hedgehog” defects 0; the loop-proliferation transition in 
certain 3D classical loop models mm ; and a columnar order¬ 
ing transition in the classical cubic dimer model (CDM) (8j- 
□1). The NCCP 1 theory involves a noncompact U(l) gauge 
theory minimally coupled to an SU(2)-symmetric complex 
field and driven through an Anderson-Higgs transition; order 
parameters for the transitions can be constructed as combina¬ 
tions of these fields. While considerable debate remains about 
the true nature of the phase transitions, and many aspects are 
not yet understood satisfactorily IE ESI, it seems clear that 
there exists at least a large range of length scales over which 
their properties are well described by the NCCP 1 model. 

Based on this viewpoint, we exploit universality to provide 
results for quantities of central importance for the transition 
in the JQ model, using numerical simulations performed on 
the CDM. The quantities in question are the scaling dimen¬ 
sions, or equivalently renormalization-group (RG) eigenval¬ 


ues y q , of operators that insert magnetic monopoles of charge 
q. Monopoles are absent in the NC CP 1 theory—according 
to an argument due to Polyakov na. they always occur at 
nonzero density in a compact U(l) gauge field but are absent 
when the field is noncompact—but provide a diagnostic for 
the transition: An oppositely charged pair of test monopoles 
is deconfined in the Coulomb (disordered) phase, but becomes 
confined in the Higgs (ordered) phase. 

In the CDM, monopoles are simply defects in the dimer 
constraint, at which dimers overlap or a site is unoccupied, 
and which carry charge in an effective gauge-theoretical de¬ 
scription mm. In the JQ model, by contrast, they corre¬ 
spond to topologically nontrivial hedgehog configurations (in 
3D space-time) of the degrees of freedom. While such de¬ 
fects generically occur at nonzero density, they are in some 
cases irrelevant at the transition, which can therefore be de¬ 
scribed by a noncompact (monopole-free) critical theory. This 
comes about because of quantum Berry phases that are asso¬ 
ciated with hedgehogs and endow monopoles with nontrivial 
transformation properties under the lattice symmetries. (They 
can therefore be associated with VBS order; see Ref. |3] and 
Section [III]) Monopoles of charge q are therefore suppressed 
at the (symmetric) critical point, unless q is an integer mul¬ 
tiple of q m determined by the rotation symmetry of the lat¬ 
tice. In order for the transition to be described by the NCCP 1 
theory, in which all monopoles are forbidden, these remain¬ 
ing monopoles must be irrelevant (i.e., y q < 0 for all nonzero 
q e q tmn Z) at the appropriate RG fixed point. 

It is therefore desirable to have a method of finding the scal¬ 
ing dimensions of charge-^ monopoles at the NC CP 1 critical 
point. Due to the topological character of hedgehogs, these 
quantities are difficult to extract directly in quantum spin mod¬ 
els. So far, only indirect evidence for their relevance or irrel¬ 
evance has been found by determining the order of the transi¬ 
tion in the JQ model on lattices of different symmetries m. 
and by studying the scaling of powers of the VBS order pa¬ 
rameter HD- The scaling dimensions can also in principle be 
calculated by generalizing the spin model to SU(/V) symmetry 
and performing a large-/V expansion, but the calculations are 
technically demanding and only results to low order in N 1 
are available mm. 

Because monopoles have a simple representation in the 
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classical CDM, and because the latter model is particularly 
amenable to Monte Carlo (MC) simulations, it provides an ef¬ 
ficient route to calculation of the quantities y q for the NCC7 J 
universality class. The goal of this work is to demonstrate that 
y q can be found using such simulations and to give explicit 
values, amounting to the first direct numerical calculation of 
these exponents. 


Outline 


In Section [IT] we introduce the cubic dimer model and de¬ 
scribe its phase structure. In Section III we review the the¬ 
ory of the Neel-VBS transition in quantum spin models, em¬ 
phasizing the effect of monopoles on the critical properties. 
Section[IV]describes the scaling of the monopole distribution 
function in the dimer model, and how this can be determined 
in numerical studies. We then describe, in Section[VJ the MC 
algorithms that we use to study the behavior of monopoles of 
charge -<7 in the dimer model. We present our numerical results 
in Section|VT|before concluding, with discussion of the impli¬ 
cations for the critical properties of quantum spin models, in 
Section En] 


the configuration energy is 

E = v 2 N 2 + V 4 IV 4 . (4) 

The continuous transition of interest occurs for vt < 0 and 
V4 > 0; in the following we choose units in which V2 = — 1 . 
(The transition is also present for V 4 = 0, but it is then found 
to show mean-field critical exponents, likely as the result of 
proximity to a multicritical point 0 .) We use a lattice with 
periodic boundary conditions and an even number L of sites 
in each direction. 


B. Phase structure and critical theory 

At high temperature T, the dimer model exhibits a Coulomb 
phase, in which occupation-number correlations are algebraic 
and test monomers are deconfined. A continuum description 
for this phase can be found by first defining an effective mag¬ 
netic field Gil El 


B p (r) = i]f 



(5) 


II. CUBIC DIMER MODEL 


where = (— is +1 on the two sublattices. The lattice 
divergence, defined by 


In this section, we briefly review the definition, phase struc¬ 
ture, and critical properties of the CDM. 


obeys 


div ? £ = ~ B ^~ 4 ) 1 . 

A* 


( 6 ) 


A. Definition 


div/fi = ii,{n(?) - 1], (7) 


Our numerical studies treat a classical statistical model of 
hard-core dimers on the links of a cubic lattice. The dimer 
occupation number d p (r) e { 0 , 1 } is defined as the number of 


dimers on the link joining the site r to its neighbor r+5 p in the 
direction p (where d fl are the unit vectors in the 3 directions). 
Apart from the defect sites (see Section II Ci, the number of 
dimers touching site r. 


n(r) = + d M (r - 5,,)], ( 1 ) 


is constrained by n(r) = 1 . 

To study the transition to an ordered state, we introduce 
interactions that count the number of nearest-neighbor parallel 
dimers. 


N 2 = zz d p {r)d p (r + <?„), ( 2 ) 

f p 

v±n 


and the number of parallel dimers around cubes of the lattice. 


^ = H d fl (r)d p (j + 6 v )d ll (r + 6 p )d p (r + 6 V + 5 p ); (3) 
f p 

v±H 

P*Pf 


and so vanishes in configurations obeying the dimer con¬ 
straint. Defects in the constraint act as charges, or magnetic 
monopoles, under this discrete Gauss law, with sign depend¬ 
ing on the sublattice. 

The partition function Z. for the dimer model, subject to the 
constraint, can be written in terms of B as 

Z = Yj^ E,T ’ ( 8 ) 

div ? 5=0 

where the energy £ of a configuration is expressed in terms of 
B. 

In a long-wavelength description, B p (r) is replaced by a 
continuum vector field B(r) with vanishing divergence, V ■ B = 
0. The effective action density in the Coulomb phase is 

mm 


Sgauge = \k\B\ 1 = ] -K\VxA\ 2 , (9) 

where /? = V x /f, and higher-order corrections, redundant 
under the RG, have been omitted. 

At a critical temperature 7 C (r’ 4 .) the system orders into a 
columnar dimer crystal, which breaks both lattice translation 
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and rotation symmetries. A suitable order parameter is the 
“magnetization” m, defined by 

m ii = ^(-l)' v ^(r). (10) 

In this ordered phase, the (connected) correlations are short- 
ranged and test monopoles are subject to confinement. 

The simultaneous appearance of order and confinement can 
be described by adding vector matter fields ifi to the gauge 
theory, so that the action density becomes 


^critical — *-*gauge ^matter 
^matter = s\(p\ 2 + u{\ip\ 2 f 
^matter-gauge |(V-iA)^| 2 , 




( 11 ) 


with pure gauge part given by Eq. ©■ The matter field tp is 
minimally coupled to A and so, being dual to the magnetic 
monopoles, can be viewed as carrying electric charge. The 
transition occurs when s is tuned through its critical value s c . 
For s < s c , <p condenses and A acquires an effective mass term 
by the Anderson-Higgs mechanism. This in turn eliminates 
the algebraic correlations and confines the monomers (mag¬ 
netic charges) through the Meissner effect. 

In the case of the dimer model, the matter field ip has 2 
components mm, and the critical action is symmetric under 
SU(2). The vector structure of the order parameter is encoded 
in the SU(2) vector structure by m = <p' dp, where cry, is a 
Pauli matrix. Replacing *S mat ter by a fixed-length constraint on 
p, one arrives at the standard form for the NC CP 1 theory. 


C. Defects at the critical point 


This quantity is related to the effective interaction V q n be¬ 
tween the pair by G q (r) = exp (-y| ff (r)/T). 

For T < T c , test charges are confined: G q (r) decreases ex¬ 
ponentially for large separation r (and any q), corresponding 
to a confining interaction, V| ff oc |r|. For T > T c , V| ff (r) has a 
finite limit, G q (r) approaches a nonzero constant, and charges 
are deconfined. 

The behavior of monopoles at the critical point can be un¬ 
derstood by considering a real-space renormalization proce¬ 
dure that preserves the discrete nature of the charges. Fol¬ 
lowing a similar logic to the case of unit-charge monopoles 
Emm, one can identify a set of scaling dimensions y q cor¬ 
responding to monopoles of charge ±q. This implies that, at 
the critical point, the monopole distribution function takes the 
form 


G q (R)~\RT 2x \ (14) 

for large separation R, where x q = d - y q . 

It should be noted that multiple defects in close proximity 
are expected to act, with regard to long-wavelength properties, 
exactly as a single defect with the same total charge. 


III. DECONFINED CRITICAL POINT IN QUANTUM 
ANTIFERROMAGNETS 


In this section, we briefly review the critical theory for the 
Neel-VBS transition in spin-j quantum antiferromagnets |3j- 
Consider the Hamiltonian 

<H = jY J 3r-3r'+ f H Q , (15) 

<rS) 


A defect in the dimer model is a site r that is touched ei¬ 
ther by no dimers or by more than one, and so has n(r) 4 1. 
In terms of the effective magnetic field, such a defect has 
div,- B 4 0 , and so can be viewed as a magnetic monopole. 
According to Eq. (jTJ, empty sites and those where two dimers 
intersect both have unit charge, with sign that alternates on the 
two sublattices. When defects occur at finite density, they de¬ 
stroy the topological order, and lead to a conventional phase 
transition between ordered and disordered phases m. 

By contrast, a test pair in an otherwise defect-free config¬ 
uration can be used as a diagnostic of the phase transition, 
using the concept of confinement. This can be given a precise 
definition in the dimer model by considering first the partition 
function evaluated in the presence of a set of charges at fixed 
positions, 

Z[QA = JV £/r - (12 > 

\B^} 

di WjtB=Qi> 


The distribution function for a test pair of monopoles of 
charge ±<7 at r ± is then defined by the ratio of the partition 
function in their presence to that in their absence. 


G q (r+ - r_) = 


Z\rj. ] 

z 


(13) 


where J > 0, the sum is over neighboring pairs of sites r and 
r' of a 2D lattice, and S, is a quantum-mechanical spin. The 
term 'Hq contains competing interactions, which can drive a 
zero-temperature transition into a VBS state, in which there is 
no magnetic order; an example is the JQ model introduced by 
Sandvik 0. 

When the first term in Eq. © dominates, the ground state 
is a Neel antiferromagnet, with a two-sublattice collinear or¬ 
dering and order parameter .) 4 0 , where e r = +1 

on the two sublattices. When 'Hq dominates, the ground 
state is instead a VBS, breaking the discrete lattice symme¬ 
tries while preserving SU(2) symmetry. An order parameter 
for this phase is the (complex) expectation value of the VBS 
operator i// V hs 0 ; different discrete ordering patterns are dis¬ 
tinguished by its phase. 

Studies of the square-lattice JQ model using quantum 
Monte Carlo mm give evidence for the claim of a direct 
continuous transition between Neel and VBS phases. This 
transition is believed 0 to be described by the same NC CP ] 
critical theory introduced in Section[II] An argument for these 
claims is sketched below; readers are referred to Ref. [3J for 
details. 

A path-integral representation of a spin- \ antiferromagnet 
can be expressed in terms of a unit vector h r ~ e r S It 
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contains, as well as terms corresponding to those in "H, a 
Berry phase l25l depending on the topology of the (periodic) 
path n r (r) in imaginary time r; its effects will be addressed 
shortly. The path integral is more conveniently written in 
terms of spinons z r , defined by h r = z' r <fz r with the constraint 
Z. 'i Z.r = 1 required for normalization of the unit vector h r . This 
definition involves a gauge redundancy under phase rotation 
Z r —» z r e I<s at any site r, and so a long-wavelength theory also 
involves a compact U(l) gauge field a. 

The Neel phase of the spin model, in which (h) + 0, 
is represented in these terms by an Anderson-Higgs phase, 
schematically expressed as a condensate of spinons, ( z r ) + 0. 
The Coulomb (deconfined) phase of the NCCP 1 theory would 
correspond to a U(l) quantum spin liquid with deconfined 
spinons. In fact, no such spin-liquid phase exists in the spin 
model, because the compact nature of the gauge field a al¬ 
lows for the existence of magnetic monopoles lfl6l . In the 
absence of a spinon condensate (i.e., in a non-Higgs phase), 
monopoles are deconfined and replace the Coulomb phase by 
a “monopole plasma”, with short-range correlations and no 
topological order. This phase, which is equivalent to the cu¬ 
bic dimer model in the disordered phase at nonzero fugacity, 
describes the VBS phase of the quantum antiferromagnet (de¬ 
spite the superficial similarity of the VBS to the dimer crystal). 

To understand the latter identification, as well as the critical 
properties, it is necessary to consider the effect of the Berry 
phases in the path integral. As shown by Haldane l26l . the 
only important contributions to the global Berry phase are as¬ 
sociated with “hedgehog” singularities of h in space-time or, 
equivalently, with monopoles. One can show ||3] that an oper¬ 
ator that inserts monopoles into the partition function has the 
same transformation properties under spatial symmetries as 
the VBS operator i/'vbs- It follows that the long-distance limit 
of the (q — 1) monopole distribution function vanishes in a 
spatially symmetric state, and conversely that a state with a 
nonzero limit must have broken spatial symmetry. These two 
statements imply, respectively, the following: (1) at a contin¬ 
uous transition into a VBS (at which point ((// V bs) = 0), the 
monopole distribution function has a zero long-distance limit; 
and (2) in the “monopole plasma” on the non-Neel side of the 
transition, there must be VBS order. 

One therefore has a scenario where single monopoles, 
though proliferating in the VBS phase, are absent at the tran¬ 
sition. The same logic can be applied to monopoles of any 
charge q for which the insertion operator transforms nontriv- 
ially under the lattice symmetries. Considering the symme¬ 
tries of the VBS operator, one finds that this argument applies 
for all q < q m j n , where q m ; n = 4 for the square lattice, q m j n = 3 
for honeycomb, and <y mm = 2 for rectangular (twofold rotation 
symmetry). Monopoles of charge q m i„ (and integer multiples 
thereof) are allowed at the critical point and will proliferate, 
leading to confinement of spinons, if the corresponding oper¬ 
ators are relevant in the RG sense. If no such monopoles are 
relevant, the effective theory describing the transition will be 
noncompact, and given by the NC CP 1 critical theory intro¬ 
duced in Section lU 

Available numerical evidence suggests that the JQ model 
has a continuous transition on the square lattice but a first- 


order transition in the case of rectangular symmetry, implying 
that monopoles of charge q — 2 are relevant at the NCCP 1 
fixed point, while those with charge <7 = 4 are irrelevant mi- 
For the honeycomb lattice, the picture is less clear l20l [321133] , 
and we therefore focus on the corresponding value of <7 = 3. 


IV. SCALING THEORY FOR MONOPOLES 


A. Monopole distribution function 


For a system of finite (linear) size L, the scaling of the 
monopole distribution function G q , given in Eq. (14 1 , can be 
extended to ! 


G q (R.L) « L~ lx '<T q (RIL) for \R\ » a, (16) 

where T q is a universal function (for each q) and a is the lattice 
scale (set to unity elsewhere but made explicit in this section). 

On the other hand, for /?| much smaller than L, whether 
larger than a or not, one expects 


G q (R, L) « g q (R) for \R\ « L , (17) 

independent of L. This follows from the definition of G q in 
terms of a charge-neutral, and hence topologically trivial, en¬ 
semble, whose free energy relative to the zero-monopole en¬ 
semble is independent of L. 

Combining Eqs. © and ( [T7] > in the intermediate regime 
a <s |/?| <sc L gives 

r,(p) ~ \,T lx « for |p| « 1 . (18) 

Together with Eq. {T6}, this simply states that the monopole 
distribution function is a power law in |/?| in the thermody¬ 
namic limit. 

The Markov-chain MC method can give the ratio of par¬ 
tition functions for two charge configurations {Q?}, provided 
that one can transition between the two with rates respecting 
detailed balance. We are therefore able to determine the ratio 


Q q (M'\G) 


G q (R, L) 
G q {R',L) 


(19) 


by implementing a MC update that allows charge -<7 


monopoles to move, described in Sections V B and V C 


A straightforward way to extract x q from MC results for Q q 
is to choose fixed R/L and R'/a, both of order unity fl4]|24i . 
Their ratio then scales as 




L- 2 ^T q (it/L) 

gq(R') 


( 20 ) 


allowing x q to be determined through finite-size scaling. Plots 
of L 2x iQ q (R,R'\L) versus p = R/L for different L (with fixed 
R') should collapse onto a single curve (for each < 7 ), propor¬ 
tional to the universal function T q . 
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An alternative method is to find Q{bR , R\ L ) for b greater 
than but of order unity. For intermediate separations, such 
that a <s |/?| <s L/b, one finds, using Eqs. (16 1 and (181, 


§(bR, R' L) 


r q (btf/L) - 

—— - » b q 

r q (R/L) 


( 21 ) 



Charge 2 
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Charge 2 


Charge 3 


Charge 3 


This form may be more reliable in cases where finite-size scal¬ 
ing is problematic, as reported in Ref. [7J 


FIG. 1. Illustration of a local charge q , consisting of a site with q + 1 
intersecting dimers. Charge q = 2 and 3 occur in the configurations 
shown above and their rotations. 


B. Additional interactions between monopoles 


V. NUMERICAL METHODS 


In the cases of most interest, where y q - d - x q is small, 
the ratio in Eq. ( [20] ) decreases rapidly with L and the relative 
statistical error increases. To improve convergence, one can 
add an explicit interaction between monopoles, by defining 

ZdQf, ( b(R)] = e ln ^ e _£dimer/T , (22) 

di\? B=Q? 


where the summation Y!?? excludes r = f and each pair of 
sites is counted only once. The function <b/ should be sym¬ 
metric and should depend on L such that it respects the pe¬ 
riodic boundaries, but is otherwise arbitrary. Ratios of parti¬ 
tion functions can be calculated by including the interaction 
<!>/ within the MC weights. (Note that <!>/ > 1 increases the 
Boltzmann weight of a configuration with oppositely charged 
monopoles.) 

The quantity Q q {R,R'\ L\ <!>/ ) can be defined by analogy to 
Eq. ([19]), and is given by 


g q aU'\L- <s> L ) = 


' <!>,(/?) ' 


S q (R, R' \ L) ■ 


(23) 


If one fixes p — R/L and R' /a of order unity and chooses 
a function <!>/ such that <\>/(pL) ~ cL 2e ^ ql (for fixed p) and 
0/ (R') ~ c'\R'\ 2e/q2 (independent of L) for some 6, then one 
finds 


Q q {R, R’\L\ ®/,) 


L 2 ^ e - d+ y^r q (R/ L ) i c w 

gq {H')\R’\ 2d ' c '’ 


(24) 


Scaling with L then gives access to the RG eigenvalue y q . The 
intermediate functional form of <!>/ should be chosen so that 
the distribution is reasonably flat, to maintain ergodicity and 
efficiency of the algorithm. 

A possible choice of additional interaction is given, for suf¬ 
ficiently large 0 , by 


®l(R) = 


V-l/9 2 


2 |* 

,fepbc(L) 



(25) 


where pbc(L) is the set of all vectors 7 such that translation by 
T is equivalent to the identity transformation under the bound¬ 
ary conditions. For |/?| <s L, the term with T — 0 dominates 
the sum and <&l(R) ~ \R\ 2e ^ r , while for R = pL with \p\ fixed 
and of order unity, <\>i(R) oc L 2e/q ~. 


This section describes the numerical methods used to cal¬ 
culate the distribution function G q for a pair of monopoles of 
charge ±q at the critical point (for q — 2 or 3). The allowed 
configurations in the constrained cubic dimer model are those 
where every site has exactly one dimer connected to it. In or¬ 
der to study the correlation function of a pair of charges, we 
consider a new set of configurations, but now with the dimer 
constraints violated at any two sites r\ and A, located on oppo¬ 
site sublattices A and B. As illustrated in Fig. |T| q + 1 dimers 
overlap at these sites, creating a charge of ±q , according to 
Eq. 0. The configuration space 

C q = \Jc q (fu? 2 ) ( 26 ) 

r { eA 

eB 


thus contains all possible dimer configurations, with all pos¬ 
sible locations r 12 of the charges. We use MC methods to 
sample with Boltzmann weight oc exp ( -E/T) over C q , where 
the energy £ of a configuration is given by Eq. 0. Given such 
samples, the pair distribution function is calculated using 


Gq(n , A) or 


number of samples in C q (r\ , A) 
number of samples in C q 


(27) 


In order to sample from the full configuration space, we 
employ two different MC update steps, both of which satisfy 
detailed balance. The first update process, T\, changes the 
configuration of the dimers but preserves the locations of the 
charges. The second update process, 7A, moves one of the 
two charges to a nearby point on the same sublattice. Details 
of the updates are given in Sections V A - V C At each MC 
step, denoted T , either T\ or 7A is chosen with equal proba¬ 
bility. Thermalization from an initial ordered state is assumed 
to have been achieved after attempting 1000 x L 3 dimer flips 
using T| and Ti. After thermalization, the average number of 
update steps Np required for L 3 dimer flips is estimated by 
averaging over several update steps. Having estimated Np, 
samples are then taken once every Np update steps T. Mul¬ 
tiple runs give independent estimates of G q allowing error es¬ 
timation using the jack-knife method. 


A. Updates of the background dimer configurations 

Here we describe the MC steps T) which sample within the 
subset C q (r\, A) with the charges located at fixed sites A,2- 
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FIG. 2. Initial step of the update process for background dimers. The 
link-monomer and the monomer at the starting node Ty are indicated 
by red/blue dots. The arrows in the figure indicate the direction of 
the link monomer. 



FIG. 4. An example of a termination step for the background dimer 
update. When the link monomer returns to the starting point To of the 
update loop, it can annihilate the monomer on the starting site and 
result in a new configuration in C q {f i, T 2 ). 



FIG. 3. Examples of individual hopping steps of the link monomer 
(blue dot). The first two rows show the scenario when the node ahead 
has no charge. The third and fourth lines show the step where the 
node ahead has a charge 2. Only hopping steps that do not change 
the local charge are allowed. These steps can be generalized to the 
case where the node has higher q. 


Any move that results in a change in the charge at r, is not 
considered, unless = r {) . In this case, the link monomer 
can also annihilate the monomer at ro, as illustrated in Fig. [f] 
completing the update step and giving a new configuration in 
C ? (n,r 2 ). 

When the link monomer hops over a site with a charge q 
(see Fig. [3j, the q + 1 dimers composing the charge are rear¬ 
ranged, but the location and charge q are unchanged. 

The probability of transitions at each individual hop of the 
link monomer, as well as of the starting and terminating steps, 
are given by 


p(k 


Wq - 8 kq min(w) 
£ vt> - min(vv) 


(28) 


where k and q represents the initial and final configurations. 
The weight of a configuration q is given by w q = exp (-E q /T), 
where E q is calculated by taking the interaction energy of a 
half dimer to be half that of a full dimer. The notations min(vv) 
and 2 w represent the minimum and the sum of the weights 
associated with all allowed final configurations q. This choice 
of probabilities ensures that detailed balance is satisfied for a 
complete step Ti. 

Note that in order to evaluate p(k —> q) only the energy dif¬ 
ferences between configurations need to be calculated. These 
energy differences can be calculated by accessing the dimer 
occupancy in a very small portion of the whole system. 


The method is based on the directed-loop algorithm |30l , and 
the specific implementation is very similar to the one used in 
Ref-GEi Here we give a summary for completeness. 

The update starts by choosing a starting node r (l , from any 
of the charge-free lattice sites, with equal probability. The 
node is accepted with a probability that depends on the local 
configuration (described below). If the step is accepted, the 
system transitions into an intermediate configuration in which 
there is a charge ± 1 (monomer) located at ro and another one 
of charge +1 on a link connected to it. Creation of two such 
charges is accompanied by addition or removal of half a dimer 
as shown in the example in Fig. [2] The link monomer has a 
‘direction’ associated with it, which initially, is set to be away 
from ro. The direction identifies one of the two nodes in that 
link as the node ‘ahead’ of the link monomer. 

In further steps, the link monomer can hop to one of the 
six links connected to the node r, ahead of it, erasing or cre¬ 
ating dimers in the process, as shown in Fig. [3] The direction 
of the link monomer after the hop is set to be away from r,. 


B. Updates for charge-2 monopole 

Here we describe the update steps T 2 used for achieving 
transitions between configurations with different locations of 
the q - 2 charges. At each step of the update one of the two 
charges can move to one of the nearest points on the same 
sublattice through local rearrangement of dimers, as shown in 

Fig-0 

In order to sample the distribution C q , the probabilities 
of the transitions need to be chosen correctly. This can be 
achieved using an acceptance-rejection method, but the differ¬ 
ent number of possible moves associated with different charge 
configurations means that the standard Metropolis probabili¬ 
ties must be modified, as follows: At the beginning of each 
step, one of the two charges is chosen at random (with equal 
probabilities). Let hq be the number of possible ways in which 
the selected charge can hop. One of these hq transitions is cho¬ 
sen at random (with equal probabilities). The selected transi¬ 
tion is accepted with probability min (1, j^p), where n\ is the 
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FIG. 5. Examples of the update step that allow a charge-2 monopole to be moved. The location of the charge is marked with red dot. Each 
possible translation of the charge involves moving exactly two of the three dimers constituting the monopole. Depending on the orientation of 
the dimers at the charge, the number of ways to hop out of the initial state can differ. 




FIG. 6. The MC procedure can be thought of as forming a graph 
with edges connecting the configurations that are related by a single 
r 2 MC step. Thus transition probabilities p a ^,t and are nonzero 
iff they are adjacent on this graph. The number of transitions out of 
(or into) a configuration a gives the degree n a . Transition probabili¬ 
ties need to be picked such that configurations are sampled with the 
desired Boltzmann probabilities. 


number of possible ways the charge can hop out of its new lo¬ 
cation, and W()j are the Boltzmann weights of the old and new 
configurations. 

To clarify the choice of probabilities, consider an initial 
configuration a with Boltzmann weight w a . Given a choice 
of one the two charges, one of the n a possible transitions is 
first chosen with probability jj -; see Fig. | 6 ] This transition is 
accepted with probability (supposing Wbn a < w a rib). The 
net effective probability of transition is thus p a -*b — jp X 
while, for the reverse process, the effective probability is 
Pb-> a = — X 1. These satisfy detailed balance w a Pa-^b = 

fib 

WbPb^a- 

Note that the number of ways no a given charge can hop 
depends on the configuration of the < 7+1 dimers connected 
to the charge. When the three dimers forming the charge (for 
q = 2 ) are coplanar, there are two ways in which the charge 


can hop to a new point. If instead the dimers are non-coplanar, 
the charge can hop in three different directions (see Fig. [5]). In 
addition, when the two charges are close to each other, the 
dimers attached to one charge can block some of the tran¬ 
sitions out of the current state, thus affecting the number of 
transitions for a given charge. 

C. Updates for charge-3 monopole 

The update steps 7~2 for transitions between configurations 
with different locations of the charge-3 monopole are simi¬ 
lar to those described for charge 2. There are zero possible 
transitions into or out of the coplanar q = 3 dimer config¬ 
uration shown in fourth panel of Fig. [T] When the dimers 
are not coplanar, (third panel of Fig. [TJ, there are two or four 
transitions possible as described here. In this configuration, 
three out of the four dimers lie on mutually perpendicular di¬ 
rections. There are two such triplets among the four dimers. 
Dimers in each such triplet lie along three adjacent edges 
defining a cube as shown in Fig. [7] If there is a dimer occu¬ 
pying an edge diagonally opposite to any of the three dimers, 
the dimers within the cube can be rearranged in two differ¬ 
ent ways as shown in the example. Each such rearrangement 
moves the charge to a new location. Since there are two dif¬ 
ferent triplets and associated cubes, there can be two or four 
transitions in total. Similar to the case of charge 2, the number 
of transitions out of a charge are again modified when the two 
charges are in close proximity. 

Note that, since there are no transitions of the type T 2 that 
can move a charge starting from a coplanar configuration, one 
would expect (from detailed balance) that there are no 7~2 tran¬ 
sitions into a coplanar configuration; this is in fact true for 
the scheme described. Ergodicity is maintained, however, be¬ 
cause a 7j update step can connect a coplanar configuration 
and a non-coplanar configuration. 
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FIG. 7. Examples of the update step that allow a charge-3 monopole 
to move. When the dimers at the charge q = 3 are in a non-coplanar 
configuration, three out of the four dimers lie on adjacent edges of 
a unit cube. There are two such triplets and two associated cubes. 
When there is a dimer at an edge diagonally opposite to any of these 
three dimers, as shown above, there are two possible transitions of 
the charge out of the configuration. Note that such transitions always 
result in a non-coplanar arrangement of dimers. 


obtained this way are identical to that obtained with full MC 
sampling in the entire region. 


2. Repulsive interaction between charges 


Another way to improve the estimation is to use an added 


repulsive potential, as explained in Section IV B Samples 
are obtained from the Boltzmann distribution where the en¬ 
ergy functional has an added repulsive interaction energy be¬ 
tween the charges of the form of Eq. ©. This interac¬ 
tion term has a sum over the repulsion from mirror charges 
in copies of the whole lattice system arising from periodic 
boundary conditions. Computationally this interaction is cal¬ 
culated by summing over 20 periodic lattices, corresponding 
to -20 L < Tj < 20 L for i = 1,2,3 in Eq. ([25). The specific 


results presented here for charge 3 were estimated with 6 — 
Estimates made with the added repulsive interaction match the 
values obtained without it within error bars. 


VI. RESULTS 

In the following subsections we describe our results for the 
monopole distribution function G q for charges q — 2 and 3. 


D. Improving convergence 

As noted in Section [IV B| the correlation function decays 
rapidly with |k| and a large number of iterations is therefore 
needed to obtain accurate estimates. A sample at a large dis¬ 
tance R ~ OiL) is obtained only once in 0(L lx 'i) steps. We 
have used two methods for improving convergence for the 
charge 3 systems. 


1. Piecewise estimation 

One way to reduce the computational time is through es¬ 
timating the correlation functions in two segments. The 
two correlation functions Q^R, R' = 1. L) and G>(R- R' = 
=R max ,L) are calculated in separate MC runs by restricting 
the samples to 

C< =C q (o<\A -ri\ < fflmax) (29) 

C> = C q Qflmax < In - ?l\ < ^max) , (30) 

respectively, where R max ^ is the maximum possible dis¬ 
tance between the charges. This is in practice achieved by 
rejecting any transition 77 that takes the system out of the 
subsets. 

The two results can be patched together at R — \R mdx to 
obtain the full function Q q . Note that there is no constraint 
on the relative locations of charge-one monomers that occur 
in the directed-loop algorithm, step T). We have checked in 
systems of sizes up to L — 30 that the correlation functions 


A. Charge 2 

Fig.[ 8 ]shows the pair correlation function Q q for charge q = 
2, rescaled by multiplying by L 2xi and plotted against R/L. 
The scaling dimension X 2 can be estimated from the plot in 
Fig. [9] showing Q^iR = pR max ,R' = 1, L) as a function of 
L for p - R/L = 0.95, using Eq. ( |20) . The calculated RG 
eigenvalues are y 2 = 1.58(2), 1.54(8), and 1.48(7) for V 4 = 
1.0, 1.5, and 10.0 respectively, in agreement within error bars. 

The dashed lines in Fig. [8] show a function oc (R/L) 2x \ us¬ 
ing the values of X 2 calculated in Fig. [9] Using Eq. ( [2l) , one 
expects the slopes of the lines to match those of the scaled 
distribution function in a region with 1 <s R <s L. We find 
slightly larger values for X 2 by this method, and an estimate 
of >’2 - 1.4 obtained by analyzing the scaling with b around 
R = 0.11 L. The narrow region in which the scaling form ap¬ 
plies, particularly for the relatively small system sizes that are 
accessible here, precludes a more precise estimate using this 
method. 


B. Charge 3 


The scaling dimension of the q = 3 charge was obtained 


using both of the methods described in Section V D Results 


using piecewise estimation, presented in Fig. 10 give an RG 
eigenvalue of >’3 = 0.31(11). The calculations were performed 
for V 4 = 1.5 and the scaling dimension was obtained using the 
monopole distribution function at a distance of R = 0.95R max . 
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FIG. 8. Monopole distribution function @ q for a pair of charges with 
q = ±2. The two figures show the data for different four-dimer inter¬ 
actions v ’4 = 1.0 (top), 1.5 (middle), and 10 (bottom). The functions 
have been rescaled with L lx i where the scaling dimension x q is re¬ 
lated to the RG eigenvalue y q obtained from Fig. [9] by x q = d - y q . 
The dashed lines have slope -2x q , and are expected to be parallel to 
the scaled correlation for 1 <K R <K L. 


FIG. 9. Log-log plot of = pRmax, R' = 1, L) for various system 
sizes at the critical temperature, with v 4 = 1.0 (top), 1.5 (middle), 
and 10 (bottom) and p = 0.95. The RG eigenvalue yi is related to the 
slopes s of the lines as s = -2x q = -2 (d - y q ). 
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FIG. 10. Top: Pair correlation function for charge q = 3, using piece- 
wise estimation. As for the case of charge 2 in Fig. [8] the function 
has been rescaled to show the scaling with L. The dashed line has a 
slope of-2jc 3 . Bottom: Log-log plot of - pL,R' = 6,L)versus 
L. The scaling dimension can be inferred directly from the slope of 
the best-fit line. 


In this case, reasonable agreement is found with Eq. ( f2T| for 
5 < R < 0.25 L. 


In order to obtain a better estimate of the exponent, we per¬ 
formed the calculations with an added repulsive interaction 
for V 4 = 1.5 and v 4 = 10. Fig. 11 shows the plot of 6 f 5 for 
this case. With the added repulsive interaction, short distance 
features appear to be amplified but the tail of the function still 
scales with the an exponent close to our previous estimate. 
The RG eigenvalues obtained were y 3 = 0.28(4) and 0.20(3) 


for V 4 = 1.5 and 10 respectively, as shown in Fig. 12 


VII. CONCLUSIONS 

We have demonstrated a method for calculating scaling di¬ 
mensions of monomers with effective charge q in a classical 
dimer model. We have applied it to the columnar ordering 
transition of dimers on the cubic lattice and calculated the RG 


FIG. 11. Pair correlation function (R = pL,R' = 1, L) calculated 
with an additional repulsive interaction between the charges of the 
form given in Eq. |25| with 6 = |. 


scaling eigenvalues y q for monomers of charge q - 2 and 3. 
For q - 2, we find the values yi = 1.58(2), 1.54(8), and 
1.48(7) for v ' 4 = 1.0, 1.5, and 10.0 respectively. For q - 3, 
we find V 3 = 0.28(4) and 0.20(3) for V 4 = 1.5 and 10.0 respec¬ 
tively. For comparison, the scaling eigenvalue for monomers 
of charge q = 1 was found as yq = 2.421(8) for v 4 = 1.0 in 
Ref. [14] (where it is denoted y z ), using the standard directed- 
loop algorithm. 

The slight drifts in the scaling dimension with V 4 are con¬ 
sistent with previous results for this transition IfDl . and have 
been attributed to corrections to scaling resulting from a 
nearby tricritical point 02 - The values for the largest V4, fur¬ 
thest from the tricritical point, should therefore be considered 
most reliable. Recent work 0 has indicated that this uni¬ 
versality class exhibits unusual finite-size effects, providing 
an alternative explanation for this behavior. As noted in Sec¬ 
tion VIA we indeed observe small discrepancies between re¬ 


sults obtained with and without assuming standard finite-size 
scaling, but considerably larger system sizes are required to 
clarify this point. 


There is now quite strong theoretical and numerical sup¬ 
port for the claim that this transition is in the NCCP 1 univer¬ 
sality class, and so other transitions in the same class should 
have identical exponents. Of most current interest are quan¬ 
tum phase transitions in 2D S - 4 antiferromagnets 0 , such 
as the JQ model a , between Neel and VBS phases. The scal¬ 
ing dimensions of monopoles of charge q , and, in particular, 
whether they are relevant or irrelevant, are of crucial impor¬ 
tance for determining the fate of such transitions. For the rect¬ 
angular (i.e., with twofold rotational symmetry) and square 
lattices, the situation is relatively clear: The nature of the tran¬ 
sition in the rectangular-lattice JQ model is determined by the 
sign of V 2 - The latter is clearly positive, and so the transition 
is certainly not described by the KCCI 11 class. For the square- 
lattice JQ model, the nature of the transition depends on the 
sign of y 4 . Previous results of quantum MC simulations on 
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FIG. 12. Log-log plot of ff 3 (R = pL,R' = 1, L) versus L, for v 4 = 1.5 (left) and 10.0 (right) and p = 0.95, calculated with an added repulsive 
potential. The continuous lines show an error-weighted least-square fit. The dotted and dashed lines correspond to the cases where y 3 would 
have been 0 and 0.5 respectively. The RG eigenvalue y 3 is related to the slope j- by s = -2(d - 9 — y 3 ), where 9 = f parametrizes the additional 
interaction. 


this model ED generally suggest that the transition is indeed 
continuous. 

For this reason, we have concentrated here on y 3 , which 
is applicable to the honeycomb lattice, where the sites 
have threefold rotation symmetry and the smallest allowed 
monopole charge is q m j n = 3. The positive value of y 3 indi¬ 
cates that the Neel-VBS transition on the honeycomb lattice 
should not be in this universality class, and would most likely 
be driven first order. Quantum MC results for such systems 
l20i[32lf33l have seen no clear evidence of a first-order transi¬ 
tion, although a trend in this direction with increasing system 
size has been suggested GO). Our results indicate that the 
true nature of this transition is indeed first order, but the small 
value of y 3 is consistent with critical behavior, described by 
the NC CP 1 universality class, over a range of length scales. 

Even allowing for the large and unconventional finite-size 
effects in this system {7][l3|. it seems unlikely that y 3 would 
be more than very weakly irrelevant, which would mean that 
three-fold anisotropy should be visible over moderate length 
scales in the JQ model on the honeycomb lattice. Pujari et al. 
Il33l have recently reported evidence of such anisotropy at the 
critical point. 

While the observed order of the transition in the JQ model 
provides evidence for the sign of y q , few quantitative values 
of y q for q > 1 are available with which to compare. By 
studying the scaling of powers of the VBS order parame¬ 
ter directly in the quantum model, Harada et al. |2(JI found 
y 2 — 1.0. One can also compare with large-/V expansions 
of the CP N 1 model ED ED, which give (in our notation) 
X 2 = 0.311 N - 0.234 + 0(N ) and x 3 = 0.544A + 0(N {) ). Trun¬ 
cating the expansions at these orders and setting N = 2 gives 


= 2.612 and y 3 = 1.912, but there is of course no reason to 
expect the truncation to be reasonable in this case. 

It would also be of interest to obtain an estimate on simi¬ 
lar lines for the monopole with <7 = 4, which is of relevance 
for the JQ model on the square lattice. In this case, however, 
single-site defects, analogous to those in Fig. |T] are even more 
difficult to move around, severely limiting the accessible sys¬ 
tem sizes. One can, in principle, construct an alternative al¬ 
gorithm where each charge 4 is realized as a fusion of two 
charge-2 monopoles situated on the nearest sites of same sub¬ 
lattice. Our preliminary results for this case show significant 
direction dependence in even the largest systems we studied 
(L = 72), making it difficult to obtain any useful estimates. 
Nonetheless, the systematic decrease of y q with increasing q 
(for 1 < q < 3) and the small value of y 3 suggest that y 4 is 
likely negative, consistent with observations of a continuous 
transition of the square-lattice JQ model ED- 
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